using namespace std;
  
#include<iostream> 
#include<fstream>
#include<fstream>
#include<cstdlib>
#include<math.h>
#include<assert.h>
#include<gsl/gsl_randist.h>

#include "../flysim.h"

//********* cExperiment member functions

cExperiment::~cExperiment(){
  cerr<<"Experiment Destructor\n";
  delete [] Dish;
  delete [] PetriDishSummaryStats; 
  delete [] TargetPetriDishSummaryStats; 
  delete [] dOldParameterValues; 
  delete [] dStatsVec; 
}
 
cExperiment::cExperiment():cFlySim() { 
  try{
    mpi_rank = MPI::COMM_WORLD.Get_rank();
    mpi_ranks = MPI::COMM_WORLD.Get_size();
    if (mpi_rank==mpi_ranks-1)  {
    cerr<<"Opening "<<outputfile<<" for write\n";
    ofile.open(outputfile.data(),ios::out);
    if (!ofile.is_open()){
      cerr<<"Cannot open the outputfile "<<outputfile<<endl;
      throw "Constructor failed for cExperiment\n";
    }
	}
    // set up dynamic arrays here
    iTotalTimeOfExperiment=(iNSamples-1)*iSamplingInterval+iTimeToStartSampling;
    // these are all the petri dishes - by convention we will make 
    // petri dish #0 represent not being on a petri dish
    Dish = new cPetriDish[gciNumberOfPetriDishes+1]; 
    // these are all the flies
    cFlyCategory * genericCategory = new cFlyCategory();
    flyCategories.push_back(genericCategory);
    for(int i=0;i<gciTotalNumberOfFlies;++i){
      cFly fly(genericCategory);
      Flies.push_back(fly);
    }
    // The following didn't look correct.  It should have size ciNoStats based
    // on subsequent bound on for loops on this array
    // float PetriDishSummaryStats[gciNumberOfPetriDishes+1][1]; // stores summary statistics of interest about the petri dishes
    PetriDishSummaryStats = new float[ciNoStats];
    // float TargetPetriDishSummaryStats[gciNumberOfPetriDishes+1][1]; // This stores the target (i.e. observed) values of the statistics
    TargetPetriDishSummaryStats = new float [ciNoStats];
    dOldParameterValues = new float[gciNumberOfFlyParameters];
    dStatsVec = new float[ciNoStats];

    iTimesAFlyMoved = iTimesAMaleMoved = iTimesAFemaleMoved = 0;
    iSamplesTaken=0;
    dTime=0;
    dLastSampleTime=0;
    iMoveCounter=0;
    iNumberOfAcceptances=0;

    // We should probably think of somewhere else to do this, but it 
    // will do for now 
    // Note that the stats are all being standardised now, 
    // with mean of 0, and SD of 1;
 
    for (int i=0;i<ciNoStats;i++)
    {
      PetriDishSummaryStats[i]=0;
      TargetPetriDishSummaryStats[i]= 0;  
    } 
     // initialize the Random Number Generator
    // use the Mersenne Twister
    rg = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(rg,ranseed);
    cerr<<"RNG created of type "<<gsl_rng_name (rg)<<
    " using seed "<<ranseed<<endl;
  }catch (exception & e){
    cerr<<"An exception was caught: "<<e.what()<<endl;
    throw "Constructor failed for cExperiment\n";
  }
}

  

void cExperiment::run()
{
//cerr<<"running run()"<<endl;
	float dSumStatsMeans[ciNoStats];
	for (int j=0;j<ciNoStats;j++)
		dSumStatsMeans[j]=0;
  if (bDoABC) // we want to do an ABC analysis
  {
	if (mpi_rank==mpi_ranks-1)
    cerr<<"Proceeding with ABC\n";
    ABC();
  }else // we just want to simulate some data
  {
    cerr<<"Proceeding with simulations\n";
	  for (int i=0;i<gciNoOfSimulationsPerProposal;i++)
	  {
    InitializeEvolutionaryParameters(1); // use the pre-defined values
		  ResetAllPetriDishStatistics();
                  statsVector.clear();
		  dSquaredDistance=0;
		  iTimesAFlyMoved = 0;
		  iTimesAMaleMoved = 0;
		  iTimesAFemaleMoved=0;
		  
		  // simulate data
		  RandomSetup(); // set the flies up randomly
		  
		  //cerr<<endl<<gciNumberOfPetriDishes<<" ";
		  //for (int iii=0;iii<gciNumberOfPetriDishes+1;iii++)
		  //{
		//	  cerr<<Dish[iii].iHowManyMalesAreHere<<" "<<Dish[iii].iHowManyFemalesAreHere<<" ";
		  //}
		  
		 // Evolve((float)(iTimeToStartSampling),0,0); // for this bit we don't record summary stats. // PM 3/18/13
    Evolve((float)(iNSamples*iSamplingInterval),1,0); // for this bit we do  // PMedit   // PM 3/18/13
    // PM changed the next line because he changed the definition of OutputParameterValues()
		  UpdateSummaryStats();
		  OutputParameterValues(i,1);
		  ofile << " " << dSquaredDistance << "  ";
		  OutputSummaryStats();
		  ofile<<"     "<<log(Flies[0].expRateLeaveMale)*3<<" "<<log(Flies[0].expRateLeaveFemale)*3<<" "<<log(Flies[0].expRateJoinMale)*3<<" "<<log(Flies[0].expRateJoinFemale)*3;
		  for (int j=0;j<ciNoStats;j++)
			  dSumStatsMeans[j]+=((float)PetriDishSummaryStats[j]);

		  
//OutputSummaryStats();
      //ofile<<endl;
      //   OutputSummaryStatsAtEndOfSim(); // PM changed the way we output the summary stats so that they output in time-order -  he then commented out everything below this because we no longer need it.

    //MyList.PrintTheEntireList();
      /*
#ifdef USE_PAUL_LINKED_LIST
    cerr<<" writing the entire list the old way\n";
    MyList.WriteTheEntireList();
#else
    cerr<<"Dumping the entire list the STL way\n";
    list<cListElement>::iterator it;
    for(it=statsList.begin();it!=statsList.end();it++){
      it->ListWrite();
    }
#endif
    ofile<<"\nSamples taken= "<<iSamplesTaken;
#ifdef USE_PAUL_LINKED_LIST
    cerr<<" emptying the entire list the STL way\n";
    MyList.EmptyTheList();
#else
    cerr<<"Clearing the entire list the STL way\n";
    statsList.clear();
#endif */
	  }
	  cerr<<"\nFinal summary stat means: ";
	  for (int j=0;j<ciNoStats;j++)
		  cerr<<dSumStatsMeans[j]/gciNoOfSimulationsPerProposal<<" ";

  }
}

void cExperiment::RandomSetup()
{
	int i,jM,jF,k,iWhere;
	//float p;
	dTime=0;
	//assign sexes
	for (i=0;i<gciTotalNumberOfMales;i++)
	{
		Flies[i].cSex='M';
		Flies[i].iID=i;
	}
	for (i=gciTotalNumberOfMales;i<gciTotalNumberOfFlies;i++)
	{	  
		Flies[i].cSex='F';
		Flies[i].iID=i;
	}
	// randomly put flies on petri dishes
	for (i=0;i<gciTotalNumberOfFlies;i++)
	{
		//p=gsl_rng_uniform_pos(rg);
		//iWhere=(int)floor(p*(gciNumberOfPetriDishes+1));
		iWhere=gsl_rng_uniform_int(rg,gciNumberOfPetriDishes+1);
		//cerr<<"I am fly "<<Flies[i].iID<<". I am on petridish "<<iWhere<<"\n";
		Flies[i].WhereAmI=&Dish[iWhere]; // the pointer 'WhereAmI' is set to point to the address (location) of PetriDish i;
		if (Flies[i].cSex=='M'){
			Dish[iWhere].iHowManyMalesAreHere++;
			// cerr<< "There are "<<Dish[iWhere].iHowManyMalesAreHere<<" males here.\n";
		}
		else{
			Dish[iWhere].iHowManyFemalesAreHere++;
			//cerr<< "There are "<<Dish[iWhere].iHowManyFemalesAreHere<<" females here.\n";
		}
	}
	for (i=0;i<gciTotalNumberOfFlies;i++)
		Flies[i].dCalcRateOfLeaving();
 	
	// set the values of the variables that keep track of how many males and females there are on each dish
	for (i=0;i<gciNumberOfPetriDishes+1;i++)
	{
		jM=0;
		jF=0;
		for (k=0;k<gciTotalNumberOfFlies;k++)
			if (Flies[k].WhereAmI->iIndex==i)
			{
				if (Flies[k].cSex=='M')
					jM++;
				else
					jF++;
			}
		Dish[i].iHowManyMalesAreHere=jM;
		Dish[i].iHowManyFemalesAreHere=jF;
	}
 	// calculate the leaving rates for each fly to initialize the process
	for (i=0;i<gciTotalNumberOfFlies;i++)
		Flies[i].dCalcRateOfLeaving();
}


void cExperiment::PrintCurrent()
{
  int i;
  for (i=0;i<gciTotalNumberOfFlies;i++)
    Flies[i].Print(1);
  ofile<<"\n\n========='Outside world':";
  Dish[0].PrintDetailsOfThisDish();
  ofile<<"\n\n==========Petri dishes:";
  for (i=1;i<gciNumberOfPetriDishes+1;i++)
    Dish[i].PrintDetailsOfThisDish();
}

void cExperiment::OldNextEvent()
{
  int i;
  float dTotalRate=0,dTot=0;
  float p;
  
  CheckState(); // for debugging purposes we make sure all variables are correct for each petri dish
  
  // to find the fly that moves we need to calculate the total rate at which they move - we just did this in the calling function, so we could store that value and use it here if speed becomes an issue (the saving will be tiny tho')
  dTotalRate=CalcEventRate();
  
  // find out which fly moves
  p=gsl_rng_uniform_pos(rg);
  i=-1;
  while ((i<gciTotalNumberOfFlies-1)&&(dTot<=p))
  {
    i++; // increase the fly index
    dTot+=Flies[i].dLeaveRate/dTotalRate; // is it this fly that moves?
  }
  assert(dTot>=p); // if this isn't true something went badly wrong
  //if i is in the first half of the n of flies it's...male?
  iTimesAFlyMoved++;
  if( Flies[i].cSex == 'M')
  {
    iTimesAMaleMoved++;
  }else{
    iTimesAFemaleMoved++;
  }
  // store the petri dish it left
  //work out where the fly goes to
  if (i<0)
    std::cout<<"\nProblem";
  PackUpAndLeaveTown(i); 
}

void cExperiment::NextEvent()
{
  int i;
  float dTot = 0;
  float dTotalRate=0;
  float p;
  
  //CheckState(); // for debugging purposes we make sure all variables are correct for each petri dish
  
  // to find the fly that moves we need to calculate the total rate at which they move - we just did this in the calling function, so we could store that value and use it here if speed becomes an issue (the saving will be tiny tho')
  dTotalRate=CalcEventRate();
  
  // find out which fly moves
  p=gsl_rng_uniform_pos(rg);
//cerr<<endl<<(Flies[0].dLeaveRate+Flies[1].dLeaveRate+Flies[2].dLeaveRate+Flies[3].dLeaveRate+Flies[4].dLeaveRate+Flies[5].dLeaveRate)/dTotalRate<<"    ";
//    for (i=0;i<gciTotalNumberOfFlies;i++)
//        cerr<<log(Flies[i].expRateLeaveMale)*3<<" ";
//    for (i=0;i<gciTotalNumberOfFlies;i++)
 //       cerr<<log(Flies[i].expRateLeaveFemale)*3<<" ";
  i=-1;
  while ((i<gciTotalNumberOfFlies-1)&&(dTot<=p))
  {
    i++; // increase the fly index
    dTot+=Flies[i].dLeaveRate/dTotalRate; // is it this fly that moves?
  }
  //if (p-dTot>0) cerr<<"Asserting will fail: "<<dTot<<","<<p<<"\n";
  //assert(dTot>=p); // if this isn't true something went badly wrong
  // GKC: The below tolerance is too stringent. We can normalize the denom if inclined.
  //assert(dTot<=1.00000001);
  //if i is in the first half of the n of flies it's...male?
  iTimesAFlyMoved++;
  if( Flies[i].cSex == 'M')
  {
    iTimesAMaleMoved++;
  }else{
    iTimesAFemaleMoved++;
  }
  // store the petri dish it left
  giPetriDishLeft=Flies[i].WhereAmI->iIndex;
  //work out where the fly goes to
  if (i<0)
    std::cout<<"\nProblem";
  PackUpAndLeaveTown(i); 
  // record where it went to
  giPetriDishArrivedAt=Flies[i].WhereAmI->iIndex;
  //  ofile<<"\nFly left dish "<<giPetriDishLeft<<" and went to "<<giPetriDishArrivedAt;
}

void cExperiment::PackUpAndLeaveTown(int iWanderingFly)
{
  int iPresentHome=Flies[iWanderingFly].WhereAmI->iIndex,iFutureHome;
  if (iPresentHome==0) // the fly is moving from the outside environment, so will land on a petri-dish
    iFutureHome=Flies[iWanderingFly].LandOnAPetriDish();
  else // the fly is on a petri dish
    iFutureHome=Flies[iWanderingFly].LeaveThisPetriDish(iPresentHome);
  // update records for that fly and the two petri dishes involved
  Flies[iWanderingFly].WhereAmI=&Dish[iFutureHome];
  UpdatePetriDishes(&Flies[iWanderingFly],iPresentHome,iFutureHome);

  //cerr<<"\nFly "<<iWanderingFly<<" moves from "<<iPresentHome<<" to "<<iFutureHome <<" \n";
  //PrintSummaryOfPetriDishes();
  
}

void cExperiment::UpdatePetriDishes(cFly *Fly, int iOldHome, int iNewHome)
{
  if (Fly->cSex=='M')
  {
    Dish[iOldHome].iHowManyMalesAreHere--;
    Dish[iNewHome].iHowManyMalesAreHere++;
  }
  else
  {
    if (Fly->cSex=='F')
    {
      Dish[iOldHome].iHowManyFemalesAreHere--;
      Dish[iNewHome].iHowManyFemalesAreHere++;
    }
    else // we have a sex-less fly
    {
      throw "\nWe have a sex-less fly. This can't be good for us or the fly. Exit.";
    }
  }
  // update the leaving rate for flies on these petri dishes
  for (int i=0;i<gciTotalNumberOfFlies;i++)
    if ((Flies[i].WhereAmI->iIndex==iOldHome)||(Flies[i].WhereAmI->iIndex==iNewHome))
      Flies[i].dCalcRateOfLeaving();
  
}

void cExperiment::PrintSummaryOfPetriDishes()
{
  //ofile<<endl;
    ofile<<"   ";
  for (int i=0;i<gciNumberOfPetriDishes+1;i++)
  {
    assert((Dish[i].iHowManyMalesAreHere+Dish[i].iHowManyFemalesAreHere>-1)&&(Dish[i].iHowManyMalesAreHere+Dish[i].iHowManyFemalesAreHere<gciTotalNumberOfFlies+1)); // just to check nothing has gone wrong
    ofile<<Dish[i].iHowManyMalesAreHere+Dish[i].iHowManyFemalesAreHere<<" ";
  }
}

void cExperiment::UpdateList(float ddTime)// PM 3/18/13 - this argument was never used, but had the same name as a global variable
{
    int i;
	
    float NumberOfBoys = Dish[1].iHowManyMalesAreHere + Dish[2].iHowManyMalesAreHere + Dish[3].iHowManyMalesAreHere + Dish[4].iHowManyMalesAreHere;
    float NumberOfGirls = Dish[1].iHowManyFemalesAreHere + Dish[2].iHowManyFemalesAreHere + Dish[3].iHowManyFemalesAreHere + Dish[4].iHowManyFemalesAreHere;
    float dAvBoys = NumberOfBoys/4;
    //  dAvBoys = (dAvBoys-dMaleAveAve)/dMaleAveSD; //standardizing the averages for boys...  // PPMM
    float dAvGirls = NumberOfGirls/4;   // PPMM
    //  dAvGirls = (dAvGirls-dFemaleAveAve)/dFemaleAveSD; //... and girls
    float dVarBoys = ((Dish[1].iHowManyMalesAreHere-dAvBoys)*(Dish[1].iHowManyMalesAreHere-dAvBoys) + (Dish[2].iHowManyMalesAreHere-dAvBoys)*(Dish[2].iHowManyMalesAreHere-dAvBoys) + (Dish[3].iHowManyMalesAreHere-dAvBoys)*(Dish[3].iHowManyMalesAreHere-dAvBoys) + (Dish[4].iHowManyMalesAreHere-dAvBoys)*(Dish[4].iHowManyMalesAreHere-dAvBoys))/4;
    dVarBoys = (log(dVarBoys + 0.1)-dMaleVarAve)/dMaleVarSD;//standardizing the ln transformed (+0.1) variances for boys
    float dVarGirls = ((Dish[1].iHowManyFemalesAreHere-dAvGirls)*(Dish[1].iHowManyFemalesAreHere-dAvGirls)+ (Dish[2].iHowManyFemalesAreHere-dAvGirls)*(Dish[2].iHowManyFemalesAreHere-dAvGirls) + (Dish[3].iHowManyFemalesAreHere-dAvGirls)*(Dish[3].iHowManyFemalesAreHere-dAvGirls) + (Dish[4].iHowManyFemalesAreHere-dAvGirls)*(Dish[4].iHowManyFemalesAreHere-dAvGirls) )/4;
    dVarGirls = (log(dVarGirls + 0.1)-dFemaleVarAve)/dFemaleVarSD; //... and girls
    float SumProductMalesFemales = (Dish[1].iHowManyMalesAreHere - dAvBoys)*(Dish[1].iHowManyFemalesAreHere - dAvGirls)+(Dish[2].iHowManyMalesAreHere - dAvBoys)*(Dish[2].iHowManyFemalesAreHere - dAvGirls)+(Dish[3].iHowManyMalesAreHere - dAvBoys)*(Dish[3].iHowManyFemalesAreHere - dAvGirls)+(Dish[4].iHowManyMalesAreHere - dAvBoys)*(Dish[4].iHowManyFemalesAreHere - dAvGirls);
    float dCovBoysGirls = SumProductMalesFemales/4;
    
    dAvBoys = (dAvBoys-dMaleAveAve)/dMaleAveSD; //standardizing the averages for boys...  // PPMM
    dAvGirls = (dAvGirls-dFemaleAveAve)/dFemaleAveSD; //... and girls    //PPMM
    
	if (dCovBoysGirls < -dMFCovCorrection) //avoid all the weirdness with very low covariances setting a floor
        dCovBoysGirls = -8.0;
    else
        //dCovBoysGirls = (log(dCovBoysGirls+dMFCovCorrection)-dMFCovarAve)/dMFCovarAve; //otherwise, a ~normal transformation
        dCovBoysGirls = (log(dCovBoysGirls+dMFCovCorrection)-dMFCovarAve)/dMFCovarSD; //otherwise, a ~normal transformation - PM replaced this line on 1/7/13
    
    float dMalesPerMale;
    if (NumberOfBoys==0)
    {
        dMalesPerMale=0;
    }
    else
    {
        dMalesPerMale = (Dish[1].iHowManyMalesAreHere*(Dish[1].iHowManyMalesAreHere-1) + Dish[2].iHowManyMalesAreHere*(Dish[2].iHowManyMalesAreHere-1) + Dish[3].iHowManyMalesAreHere*(Dish[3].iHowManyMalesAreHere-1) + Dish[4].iHowManyMalesAreHere*(Dish[4].iHowManyMalesAreHere-1))/NumberOfBoys;
        dMalesPerMale=(log(dMalesPerMale+0.1)-dMpMAve)/dMpMSD;
    }
    
    float dFemalesPerMale;
    if (NumberOfBoys==0)
    {
        dFemalesPerMale =0;
    }
    else
    {
        dFemalesPerMale = (Dish[1].iHowManyMalesAreHere*Dish[1].iHowManyFemalesAreHere + Dish[2].iHowManyMalesAreHere*Dish[2].iHowManyFemalesAreHere + Dish[3].iHowManyMalesAreHere*Dish[3].iHowManyFemalesAreHere + Dish[4].iHowManyMalesAreHere*Dish[4].iHowManyFemalesAreHere)/NumberOfBoys;
        dFemalesPerMale=(log(dFemalesPerMale+0.1)-dFpMAve)/dFpMSD;
    }
    
    float dMalesPerFemale;
    if (NumberOfGirls==0)
    {
        dMalesPerFemale=0;
    }
    else
    {
        dMalesPerFemale= (Dish[1].iHowManyMalesAreHere*Dish[1].iHowManyFemalesAreHere + Dish[2].iHowManyMalesAreHere*Dish[2].iHowManyFemalesAreHere + Dish[3].iHowManyMalesAreHere*Dish[3].iHowManyFemalesAreHere + Dish[4].iHowManyMalesAreHere*Dish[4].iHowManyFemalesAreHere)/NumberOfGirls;
        dMalesPerFemale=(log(dMalesPerFemale+0.1)-dMpFAve)/dMpFSD; 
    }
    
    float dFemalesPerFemale;
    if (NumberOfGirls==0)
    {
        dFemalesPerFemale=0;
    }
    else
    {
        dFemalesPerFemale= (Dish[1].iHowManyFemalesAreHere*(Dish[1].iHowManyFemalesAreHere-1) + Dish[2].iHowManyFemalesAreHere*(Dish[2].iHowManyFemalesAreHere-1) + Dish[3].iHowManyFemalesAreHere*(Dish[3].iHowManyFemalesAreHere-1) + Dish[4].iHowManyFemalesAreHere*(Dish[4].iHowManyFemalesAreHere-1))/NumberOfGirls;
        dFemalesPerFemale=(log(dFemalesPerFemale+0.1)-dFpFAve)/dFpFSD;
    } 
    
    
    //float dVarLnAvBoys = 0;
    //float dVarLnAvGirls = 0;
    
    //float dStatsVec[ciNoStats] = {dAvBoys, dAvGirls, dVarBoys, dVarGirls, dCovBoysGirls, iTimesAFlyMoved,iTimesAMaleMoved,iTimesAFemaleMoved};
    // PMM
	assert(ciNoStats==10); // The code below kind of assumes there can only be 6 different statistics.  We should change below
    int j=0;
    dStatsVec[j++] = dAvBoys*dAvBoysWeight;
    dStatsVec[j++] = dAvGirls*dAvGirlsWeight;
    dStatsVec[j++] = dVarBoys*dVarBoysWeight;
    dStatsVec[j++] = dVarGirls*dVarGirlsWeight;
    dStatsVec[j++] = dCovBoysGirls*dCovBoysGirlsWeight;
	dStatsVec[j++] = dMalesPerMale*dMalesPerMaleWeight;
	dStatsVec[j++] = dFemalesPerMale*dFemalesPerMaleWeight;
	dStatsVec[j++] = dMalesPerFemale*dMalesPerFemaleWeight;
	dStatsVec[j++] = dFemalesPerFemale*dFemalesPerFemaleWeight;
    // PMM new stats  - the sd of first two stats
	//dStatsVec[j++] = dTime;
    //cerr<<endl<<dTime<<" VVV "<<ddTime;
	dStatsVec[j++] = ddTime;
	
    //cerr<<"Adding vector element\n";
    cListElement element;     // Paul wonders whether this is right. Does cListElement have the same type as a vector?
    for (i=0; i<ciNoStats; i++){
        element.dStats[i] = dStatsVec[i];
    }
    statsVector.push_back(element);

}

float dMAveSD = 0; // Temporary patch to allow us to calculate a couple of extra stats
float dFAveSD = 0;
// this is the version that causes the segmentation fault
/* void cExperiment::UpdateSummaryStats()
{
    // A list iterator to go backwards in the linked list
    //list<cListElement>::reverse_iterator rit;
    // A list iterator to go forwards in the linked list
    //list<cListElement>::iterator it;
    //cerr<<"Begin updateSummaryStats\n";
    //BRF_08_20: I added the machinery here to sample the output statistics at equally spaced intervals
    //accounting for the time between moves. I hope.
    int i, j, k;
    int iNoOfSamplesTakenFromLastInterval; // PM added this to keep track of taking multiple samples form the same interval
    float dStartSampleTime, dEndSampleTime, dSampleInterval, dSampleRemainder;
    //float dLastSampleTime;
    float dLastTime=0; // PM 4/2/13 part of a fix to make sure multiple samples can be taken from long intervals 
    //BRF_08_20: These are the variables I need to do the sampling.
    
    //BRF_09_Mar: Here I'm going to calculate standard deviation on male ave and female ave
    dMAveSD = 0;
    dFAveSD = 0;
    float dMAveVec[iNSamples];
    float dFAveVec[iNSamples];

    
    // here I'm gonna calculate the actual summary stats (not really 'petri dish' specific,
    // but we'll stick with the lingo - for now)
    k=0;
    
    for (i=0; i<ciNoStats; i++)
        PetriDishSummaryStats[i] = 0;
    
    dEndSampleTime = statsVector[statsVector.size()-1].dStats[ciNoStats-1];
    
    //cerr<<"\ndEndSampleTime "<<dEndSampleTime;
    
//cerr<<"looping from statsvector begin\n";
    vector<cListElement>::iterator vit = statsVector.begin();
    // Brad wants to skip the bottom element of the list somehow (BRF_08_20),
     vit++; 
    if (vit!=statsVector.end()){
        dStartSampleTime = vit->dStats[ciNoStats-1];
        //cerr<<"\ndStartSampleTime="<<dStartSampleTime;
        //iNSamples-1 intervals gives us iNSamples
        dSampleInterval=(dEndSampleTime-dStartSampleTime)/(iNSamples-1);
        dSampleRemainder=0;
    }
//cerr<<"looping from statsvector end\n";
    //cerr<<"\n***dSampleInterval="<<dSampleInterval<<endl;
    //for (int jj=0;jj<iNSamples;jj++)
     //   cerr<<dStartSampleTime+dSampleInterval*(jj)<<"  ";
    
    // new STL way
    // Brad wants to ignore the element at the top of the list and the
    // boundary condition below reflects this:
    unsigned int curIndex = 0;
    vit = statsVector.begin();
    //cerr<<"\n\nIterating through "<<statsVector.size()<<" vector elements";
    while (k<iNSamples){
    //while(curIndex<(statsVector.size()-1)){
        iNoOfSamplesTakenFromLastInterval=0;
        //if(dSampleRemainder <= 0) {
        if(dSampleRemainder <= 1e-9) {   // changed 0 to 1e-9 to avoid missing the last step do to rounding error. 
            ++k;
            iNoOfSamplesTakenFromLastInterval++;
            //cerr<<"\n*******Sample "<<k<<" taken*****";
            //cerr<<"\nNew dStats: ";
            for (j=0; j<ciNoStats; j++){
                PetriDishSummaryStats[j] = PetriDishSummaryStats[j] + (vit->dStats[j]);
                //cerr<<vit->dStats[j] << " ";
            }
            //BRF_09_Mar: Here I'm going to calculate standard deviation on male ave and female ave
            dMAveVec[(k-1)] = (vit -> dStats[0]);
            dFAveVec[(k-1)] = (vit -> dStats[1]);
            
            dSampleRemainder+=dSampleInterval;
            //cerr<<"    Remainder:"<<dSampleRemainder;
            //dLastSampleTime=vit->dStats[ciNoStats-1]+(iNoOfSamplesTakenFromLastInterval-1)*dSampleInterval;
            //cerr<<"\n***dLastSampleTime="<<dLastSampleTime;
       }
        //dLastSampleTime=vit->dStats[ciNoStats-1];
        if ((dSampleRemainder>0)&&(k<iNSamples)) // we skip this bit if we need to take more than one sample from this point or if we have just taken the last sample- PM 4/2/13
        {
            if (curIndex<statsVector.size()) // PM added the if statement to stop us looping past the end of the list 4/2/13
            {
                dLastTime=vit->dStats[ciNoStats-1]; // record the time of the last interval
                vit++;}
            // we're going from the last to the first sample taken
            //cerr<<"\nLine skipped: ";
            for (j=0; j<ciNoStats; j++){
                //PetriDishSummaryStats[j] = PetriDishSummaryStats[j] + (vit->dStats[j]);
                //cerr<<vit->dStats[j] << " ";
            }
            dSampleRemainder-= (vit->dStats[ciNoStats-1] - dLastTime );
            //cerr<<"\ndSampleRemainder="<<dSampleRemainder;
            ++curIndex;
        }
    }
//cerr<<"2 looping from statsvector end\n";
    //ofile << "\nI just grabbed the line: " ;
    // k++; // PM removed this on 4/2/13 becuase he has reworked this function.
    
    // PM 4/2/13 we shouldn't need the following loop any more. The last entry should already have been dealt with
   // while (vit!=statsVector.end()){
    //    vit++;
   // cerr<<"\nFinal dStats: ";
     //   for (j=0; j<ciNoStats; j++){
        // manually add in the last list entry
      //  if (vit!=statsVector.end()){
       //     PetriDishSummaryStats[j] = PetriDishSummaryStats[j] + (vit->dStats[j]);
       //     cerr<<vit->dStats[j] << " " ;
        //}
        //}
    //}
    
    //BRF_09_Mar: Here I'm going to calculate standard deviation on male ave and female ave
    dMAveVec[(k-1)] = (vit -> dStats[0]);
    dFAveVec[(k-1)] = (vit -> dStats[1]);
//cerr<<"3 looping from statsvector end\n";

    //ofile<<"k= "<<k;
    //ofile<<"\nAnd the output is: ";
    for (i=0; i<ciNoStats; i++)
    {
        PetriDishSummaryStats[i] = PetriDishSummaryStats[i]/k;
        //  ofile<<PetriDishSummaryStats[i]<<" ";
    }
    if (k!=iNSamples){
        ofile<<"\nWe didn't collect the correct number of samples to form the statistics. k="<<k;
    }
    //if (statsList.size()>20) throw "Checkpoint";
    
    
    //BRF_09_Mar: Here I'm going to calculate standard deviation on male ave and female ave
    for (i=0; i<k; i++){
        dMAveSD = ( dMAveSD + pow((dMAveVec[i]-(PetriDishSummaryStats[0])), 2));
        dFAveSD = ( dFAveSD + pow((dFAveVec[i]-(PetriDishSummaryStats[1])), 2));
        //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";
    }
    //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";
    dMAveSD = pow((dMAveSD/(k-1)), 0.5);
    dFAveSD = pow((dFAveSD/(k-1)), 0.5);
    //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";
    dMAveSD = pow(log(dMAveSD + 0.01), 2);
    dFAveSD = pow(log(dFAveSD+0.01), 2);
    //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";

    // ofile<<"\nPetriStats ";
    //for (int j=0; j<ciNoStats; j++)
    //    ofile<<PetriDishSummaryStats[j]<<" ";
    
//ofile<<"\nDumping summary stats\n";
//    rit=statsList.rbegin();
//    while (rit!=statsList.rend())
 //   {
  //       {
 //          for (j=0; j<ciNoStats; j++)
 //               ofile<<rit->dStats[j]<<" ";
 //           ofile<<"\n";
 //       }
 //       rit++;
 //   }

}*/

void cExperiment::UpdateSummaryStats()
{
    // A list iterator to go backwards in the linked list
    list<cListElement>::reverse_iterator rit;
    // A list iterator to go forwards in the linked list
    list<cListElement>::iterator it;
    //cerr<<"Begin updateSummaryStats\n";
    //BRF_08_20: I added the machinery here to sample the output statistics at equally spaced intervals
    //accounting for the time between moves. I hope.
    int i, k;
    float dStartSampleTime, dEndSampleTime, dSampleInterval;/*dLastSampleTime,*/ /*dSampleRemainder,*/ /*dLastSampleMAveVec, dLastSampleFAveVec;*/
    //BRF_08_20: These are the variables I need to do the sampling.
    
    //BRF_09_Mar: Here I'm going to calculate standard deviation on male ave and female ave
    dMAveSD = 0;
    dFAveSD = 0;
    float dMAveVec[iNSamples];
    float dFAveVec[iNSamples];
    
    // here I'm gonna calculate the actual summary stats (not really 'petri dish' specific,
    // but we'll stick with the lingo - for now)
    k=0;
    
    for (i=0; i<ciNoStats; i++)
        PetriDishSummaryStats[i] = 0;
    
    dEndSampleTime = statsVector[statsVector.size()-1].dStats[ciNoStats-1];
    
    // cerr<<"\n "<<dEndSampleTime;
    
    vector<cListElement>::iterator vit = statsVector.begin();
    // Brad wants to skip the bottom element of the list somehow (BRF_08_20),
    // so we'll do that as well here:
//    vit++;
    if (vit!=statsVector.end()){
        dStartSampleTime = vit->dStats[ciNoStats-1];
        //iNSamples-1 intervals gives us iNSamples
        dSampleInterval=(dEndSampleTime-dStartSampleTime)/(iNSamples-1);
        //dSampleRemainder=0;
    }
    
    // new STL way
    // Brad wants to ignore the element at the top of the list and the
    // boundary condition below reflects this:
    //unsigned int curIndex = 0;
   // vit = statsVector.begin();

    //cerr<<"Iterating through "<<statsVector.size()<<" vector elements\n";
    //while (vit!=statsVector.end()){
    //while(curIndex<(statsVector.size()-1)){ // seg. fault caused byfact that statVector.end() is not always defined for some reason

//alex method 5/15/13:

//cout<<"Interval is: "<<dSampleInterval<<" ,from "<<dStartSampleTime<<" to "<<dEndSampleTime<<" with "<<iNSamples<<"samples"<<endl;

dMAveVec[0]=vit -> dStats[0];
dFAveVec[0]=vit->dStats[1];
/*cout<<"Taking first sample at"<<dStartSampleTime<<"?="<< vit->dStats[ciNoStats-1] <<endl;
*/
//dLastSampleTime=vit->dStats[ciNoStats-1];
//dLastSampleMAveVec=vit->dStats[0];
//dLastSampleFAveVec=vit->dStats[1];

//float vit_tmp[ciNoStats];
//float sample_times[iNSamples];

for (int j=0; j<ciNoStats; j++) { // HOW ON EARTH DID I FORGTE THIS 5/24/... still
        PetriDishSummaryStats[j] = (vit->dStats[j]);
//cout<<"0round has PetriDishSummaryStats["<<j<<"]="<<PetriDishSummaryStats[j]<<endl;
}

//sample_times[0]=vit->dStats[ciNoStats-1];
vector<cListElement>::iterator last_sample=statsVector.begin();

for (k=1;k+1<iNSamples;k++) {
        for (vit=last_sample;vit!=statsVector.end();vit++)  {
                if ((dStartSampleTime+k*dSampleInterval) <=  vit->dStats[ciNoStats-1]) { // 
                        dMAveVec[k]=last_sample->dStats[0];//vit_tmp[0]; //vit->dStats[0];//vit_tmp[0];//dLastSampleMAveVec;
                        dFAveVec[k]=last_sample->dStats[1];//vit_tmp[1];//vit->dStats[1];//vit_tmp[1];//dLastSampleFAveVec;                     
                        for (int j=0; j<ciNoStats; j++){ // HOW ON EARTH DID I FORGTE THIS 5/24/... still
                                PetriDishSummaryStats[j] += (last_sample->dStats[j]);
                        //      cout<<k<<"round has PetriDishSummaryStats["<<j<<"]="<<PetriDishSummaryStats[j]<<" after adding "<< vit_tmp[j] <<endl;
                        }
                        //sample_times[k]=last_sample->dStats[ciNoStats-1];
                        break;
                }

                last_sample=vit;
        }
}
vit=statsVector.end()-1;
//cout<<"and out last sample at: "<<dEndSampleTime<<"?="<<vit->dStats[ciNoStats-1]<<endl;
dMAveVec[iNSamples-1]=vit->dStats[0];
dFAveVec[iNSamples-1]=vit->dStats[1];
//sample_times[iNSamples-1]=vit->dStats[ciNoStats-1];
for (int j=0; j<ciNoStats; j++){ // HOW ON EARTH DID I FORGTE THIS 5/24/... still
        PetriDishSummaryStats[j] += (vit->dStats[j]);
}



    for (i=0; i<ciNoStats; i++)
    {
        PetriDishSummaryStats[i] = PetriDishSummaryStats[i]/k;
        //  ofile<<PetriDishSummaryStats[i]<<" ";
    }
    
    
    //BRF_09_Mar: Here I'm going to calculate standard deviation on male ave and female ave
    for (i=0; i<k; i++){
        dMAveSD = ( dMAveSD + (dMAveVec[i]-PetriDishSummaryStats[0])*(dMAveVec[i]-PetriDishSummaryStats[0]));
        // GKC: The line below must be a typo! Shouldn't be dFAveVec?
        //dFAveSD = ( dFAveSD + (dMAveVec[i]-PetriDishSummaryStats[1])*(dMAveVec[i]-PetriDishSummaryStats[1]));
        dFAveSD = ( dFAveSD + (dFAveVec[i]-PetriDishSummaryStats[1])*(dFAveVec[i]-PetriDishSummaryStats[1]));
        //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";
    }
    //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";
    dMAveSD = pow((dMAveSD/(k-1)), 0.5);
    dFAveSD = pow((dFAveSD/(k-1)), 0.5);
    //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";
    dMAveSD = log(dMAveSD + EPS)*log(dMAveSD + EPS);
    dFAveSD = log(dFAveSD+EPS)*log(dFAveSD+EPS);
    //ofile <<"dMAveSD = " << dMAveSD << " dFAveSD = " << dFAveSD << " \n";

    
    //if (statsList.size()>20) throw "Checkpoint";
    
    // ofile<<"\nPetriStats ";
    //for (int j=0; j<ciNoStats; j++)
    //    ofile<<PetriDishSummaryStats[j]<<" ";
    
    //ofile<<"\nDumping summary stats\n";
    //    rit=statsList.rbegin();
    //    while (rit!=statsList.rend())
    //   {
    //       {
    //          for (j=0; j<ciNoStats; j++)
    //               ofile<<rit->dStats[j]<<" ";
    //           ofile<<"\n";
    //       }
    //       rit++;
    //   }
    
}





void cExperiment::OutputSummaryStats()
{  
  
  // UpdateSummaryStats(); // PM 3/18/13 they are updated every simulation anyway, so we don't need to do it agan here.
  //if (bUseBortotScheme)
    ofile<<"    "<<iMoveCounter<<" ";    
   // else
   //   ofile<<"Summary stats: (#events= "<<iMoveCounter<<")  ";
  for (int i=0;i<ciNoStats ;i++){
    ofile<<((float)PetriDishSummaryStats[i]) << " ";
  }
    ofile<<dMAveSD << " " << dFAveSD << " ";

  
  if (bDoABC == 0) {
    ofile <<iTimesAFlyMoved << " " << iTimesAMaleMoved << " " << iTimesAFemaleMoved << " " ;
  }
//  if (!bUseBortotScheme)
//    ofile << "\n";
  
}


void cExperiment::CheckState()
{
  int i,j,kM,kF;
  // check that the recorded number of males and females on each fly is correct
  for (i=0;i<gciNumberOfPetriDishes+1;i++)
  {
    kM=0;
    kF=0;
    for (j=0;j<gciTotalNumberOfFlies;j++)
    {
      if (Flies[j].WhereAmI->iIndex==i) // we found a fly on this dish
      {
        if (Flies[j].cSex=='M')
          kM++;
        else
          kF++;
      }
    }
    assert(Dish[i].iHowManyMalesAreHere==kM);
    assert(Dish[i].iHowManyFemalesAreHere==kF);
  }
  
  // check that we still have the correct number of flies
  j=0;
  for (i=0;i<gciNumberOfPetriDishes+1;i++)
    j+=(Dish[i].iHowManyMalesAreHere+Dish[i].iHowManyFemalesAreHere);
  assert(j==gciTotalNumberOfFlies);
}


void cExperiment::Evolve(int iHowLong, bool bRecordSummaryStats, bool bOutputSummaryStats)  //BRF_08_14
{
  int i;
  i=0;
    float dTT=0;
  float dEventRate,dNextEventTime;
        //float dT,dStartTime=dTime;
    while (i<iHowLong)// PM 3/18/13
  //  while (dTT<iHowLong)
  {
    //CheckState(); // just to make sure everything is ok with the record-keeping variables and rates
  
    // how long to the next event?
    dEventRate=CalcEventRate();
      //cout<<"\nRate="<<dEventRate;
    // ofile<<endl;
	 // for (int j=0;j<5;j++)
    //      ofile<<Dish[j].iHowManyMalesAreHere+Dish[j].iHowManyFemalesAreHere<<" ";
  
    //dNextEventTime=ExponentialRV(dEventRate); // times to events are exponentially distributed under this model
    dNextEventTime =  gsl_ran_exponential (rg, 1.0/dEventRate); // PM changed on 3/18/13 gsl_ran_exponential uses the mean as its parameter, not the rate
     // cerr<<"\ndNextEventTime="<<dNextEventTime<<"       dTime="<<dTime;
  //  PPP  ofile<<"\nFly configuration:";
  //    for (int iWhere=0;iWhere<gciNumberOfPetriDishes+1;iWhere++)
  //        ofile<<Dish[iWhere].iHowManyMalesAreHere+Dish[iWhere].iHowManyFemalesAreHere<<" ";
  //    ofile<<"\nLeaving Rates: ";
  //    for (int j=0;j<gciTotalNumberOfFlies;j++)
  //        ofile<<Flies[j].dLeaveRate<<" ";
  //    ofile<<"\ndNextEventTime ="<<dNextEventTime;
    // sample summary stats in the current state if appropriate
    //while (dLastSampleTime+dSamplingInterval<dTime)  BRF_08_14
    //{
      // update the last sample time
    dLastSampleTime+=dSamplingInterval;
      
      // update and sample the stats
    if (bRecordSummaryStats)
    {
        //if(iMoveCounter%iSamplingInterval==0)
        //{
        UpdateList(dTime);
        //UpdateList(dNextEventTime); // PM update 3/18/13
        //ofile<<"\nEvent list after update but before move:";
        //list<cListElement>::iterator it;
        //for(it=statsList.begin();it!=statsList.end();it++){
        //    it->ListWrite();}

        
//UpdateSummaryStats(); // PM
        iSamplesTaken++; // record that we have taken a sample
        if (bOutputSummaryStats)
        {
            if (bDoABC) // we want to do an ABC analysis
                OutputSummaryStats();  //working here
            else
            {
                if (iMoveCounter%iSamplingInterval==0)
                    OutputSummaryStats();  //working here
            }
        }
        //}
    }
       if (bOutputSummaryStats)
          PrintSummaryOfPetriDishes();
    dTime+=dNextEventTime;
      dTT+=dNextEventTime;
      //if (dNextEventTime>10000)
      //{
      //    ofile<<"\n***bOutputSummaryStats: "<<bOutputSummaryStats<<"       dNextEventTime= "<<dNextEventTime;
      //    PrintSummaryOfPetriDishes();
      //    ofile<<endl;
      //    for (int kk=0;kk<gciTotalNumberOfFlies;kk++)
      //    {
      //    ofile<<"     "<<Flies[i].dLeaveRate;
      //    }
      //}
      // do the next event
    NextEvent();
      // count it
    iMoveCounter++;
    i++;
      
 
  
  }
    
    // PM added this final update on 3/18/13
    if (bRecordSummaryStats)
    {
        UpdateList(dTime);
        //UpdateList(dNextEventTime); // PM update 3/18/13
        //ofile<<"\nEvent list after update but before move:";
        //list<cListElement>::iterator it;
        //for(it=statsList.begin();it!=statsList.end();it++){
        //    it->ListWrite();}
    }

    
    // if the next event happened before time ran out....
    //if (dTime+dNextEventTime<dStartTime+dHowLong)
    //{
    //  // update the clock
    //  dTime+=dNextEventTime;
    //  // do the next event
    //  NextEvent();
    //  // count it
    //  iMoveCounter++;
    //}
    //else // just need to update the clock
    //  dTime=dStartTime+dHowLong;
  
}

/*void cExperiment::Evolve(int iHowLong, bool bRecordSummaryStats, bool bOutputSummaryStats)  //BRF_08_14
{
	int i;
	i=0;
	float dEventRate,dNextEventTime;
	//float dT,dStartTime=dTime;
	while (i<iHowLong)
	{
		//CheckState(); // just to make sure everything is ok with the record-keeping variables and rates
		
		// how long to the next event?
		dEventRate=CalcEventRate();
		//dNextEventTime=ExponentialRV(dEventRate); // times to events are exponentially distributed under this model
		dNextEventTime =  gsl_ran_exponential (rg, 1.0/dEventRate);
		
		// sample summary stats in the current state if appropriate
		//while (dLastSampleTime+dSamplingInterval<dTime)  BRF_08_14
		//{
      // update the last sample time
		dLastSampleTime+=dSamplingInterval;
      
      // update and sample the stats
		if (bRecordSummaryStats)
		{
			//if(iMoveCounter%iSamplingInterval==0)
			//{
			UpdateList(dTime);
			//UpdateSummaryStats(); // PM
			iSamplesTaken++; // record that we have taken a sample
			if (bOutputSummaryStats)
			{
            if (bDoABC) // we want to do an ABC analysis
					OutputSummaryStats();  //working here
            else
            {
					if (iMoveCounter%iSamplingInterval==0)
						OutputSummaryStats();  //working here
            }
			}
			//}
		}
      
      if (bOutputSummaryStats)
			PrintSummaryOfPetriDishes();
		dTime+=dNextEventTime;
      // do the next event
		NextEvent();
      // count it
		iMoveCounter++;
		i++;
		
	}
	
	// if the next event happened before time ran out....
	//if (dTime+dNextEventTime<dStartTime+dHowLong)
	//{
	//  // update the clock
	//  dTime+=dNextEventTime;
	//  // do the next event
	//  NextEvent();
	//  // count it
	//  iMoveCounter++;
	//}
	//else // just need to update the clock
	//  dTime=dStartTime+dHowLong;
	
}*/


float cExperiment::CalcEventRate()
{
  //cerr<<"calcEventRate\n";
  float dRate=0;
  int i;
  // for each fly - find its rate of leaving its current location
  //  ofile<<"\nLeaving rates: ";
  for (i=0;i<gciTotalNumberOfFlies;i++)
  {
    // only re-calculate it if this fly is one a plate that has changed
//    if ((Flies[i].WhereAmI->iIndex==giPetriDishLeft)||(Flies[i].WhereAmI->iIndex==giPetriDishArrivedAt))
      Flies[i].dCalcRateOfLeaving();
 //     cerr<<endl<<i<<" xx "<<Flies[i].dLeaveRate<<" "<<Flies[i].WhereAmI->iIndex;
    dRate+=Flies[i].dLeaveRate;
      //if (gcPaul==1)
      //    ofile<<"\nFlychecked. Rate: "<<Flies[i].dLeaveRate;
      //ofile<<Flies[i].WhereAmI->iIndex<<" "<<Flies[i].dLeaveRate<<"   ";
    //cerr<<"Rate so far: "<<dRate<<endl;
  }
  return dRate;
}

// PM added this function
void cExperiment::ChangeTolerance()
{
  float p,q;
  p=gsl_rng_uniform_pos(rg);
  q=gsl_rng_uniform_pos(rg);
  if (p<cdProbToleranceMovesDown) // decrease it
  {
    dABCTolerance-=q*gdMaxToleranceChange;
    if (dABCTolerance<gdMinTolerance)
      dABCTolerance=gdMinTolerance-(dABCTolerance-gdMaxTolerance);
  }
  else // increase it
  {
    dABCTolerance+=q*gdMaxToleranceChange;
    if (dABCTolerance>gdMaxTolerance)
      dABCTolerance=gdMaxTolerance-(dABCTolerance-gdMaxTolerance);
  }
}

// PM changed this function
void cExperiment::ProposeNewState()
{
  int i;
  float p,q,dNewParam;
  float dUpperBound=20,dChangeAmount;
  
  // now we might change the threshold
  p=gsl_rng_uniform_pos(rg);
  bToleranceChanged=0;
  if((p<cdProbChangeTolerance)&&(bUseBortotScheme))
  {
    //change the threshold - for now we use something that tends to go downwards.
    gdOldBortotTolerance=dABCTolerance;
    this->ChangeTolerance();
    bToleranceChanged=1;

  }
  // currently we assume all flies have the same parameters - we need to change this function if not
  // First we store the old parameter values for fly movement so that we can return to them if an ABC move is rejected
  for (i=0;i<gciTotalNumberOfFlies;i++)
  {
    dOldParameterValues[0]=Flies[i].RateLeaveMale;
    dOldParameterValues[1]=Flies[i].RateJoinMale;
    dOldParameterValues[2]=Flies[i].RateLeaveFemale;
    dOldParameterValues[3]=Flies[i].RateJoinFemale;
    dOldParameterValues[4]=Flies[i].getMaleAttractionToMales(); 
    dOldParameterValues[5]=Flies[i].getMaleAttractionToFemales(); 
    dOldParameterValues[6]=Flies[i].getFemaleAttractionToMales(); 
    dOldParameterValues[7]=Flies[i].getFemaleAttractionToFemales(); 
  }
  
  // now propose new state - for now we pick a parameter at random and add a uniform [0,10] r.v. to it
  // we reflect at the boundaries of 0 and dUpperBound to make the Hasting's Ratio work out nicely
  //p=gsl_rng_uniform_pos(rg);
  // pick the parameter
  //i=(int)floor(p*(gciNumberOfFlyParameters-1)); 
  i=gsl_rng_uniform_int(rg,gciNumberOfFlyParameters-1); 
  
  // sometimes we make a small change; sometimes we make a big change  
  p=gsl_rng_uniform_pos(rg);
  q=10*gsl_rng_uniform_pos(rg)-5;
  // little av=0.05 med av=1 large av=5
  if(p<0.3)
    dChangeAmount=q/50;  
  else if (p<0.7) // make a small change
    dChangeAmount=q/10;
  else // make a big change
    dChangeAmount=q/5;
  
  
  float eparam;
  switch (i)
  {
    case 100: //this will stay 0
      dNewParam=Flies[0].RateLeaveMale+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      for (i=0;i<gciTotalNumberOfFlies;i++)
        Flies[i].RateLeaveMale=dNewParam;
      break;
    case 0:
      dNewParam=Flies[0].RateJoinMale+10+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; //BRF_08_10
      eparam = exp(dNewParam/3);
      for (i=0;i<gciTotalNumberOfFlies;i++)
      {
        Flies[i].RateJoinMale=dNewParam;
        Flies[i].expRateJoinMale=eparam;
      }
      break;
    case 1:
      dNewParam=Flies[0].RateLeaveFemale+10+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; 
      eparam = exp(dNewParam/3);
      for (i=0;i<gciTotalNumberOfFlies;i++)
      {
        Flies[i].RateLeaveFemale=dNewParam;
        Flies[i].expRateLeaveFemale=eparam;
      }
      break;
    case 2:
      dNewParam=Flies[0].RateJoinFemale+10+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; 
      eparam = exp(dNewParam/3);
      for (i=0;i<gciTotalNumberOfFlies;i++)
      {
        Flies[i].RateJoinFemale=dNewParam;
        Flies[i].expRateJoinFemale=eparam;
      }
      break;
    case 3:
      dNewParam=Flies[0].getMaleAttractionToMales()+10+dChangeAmount;
      
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; 
      //for (i=0;i<gciTotalNumberOfFlies;i++)
        //Flies[i].setMaleAttractionToMales(dNewParam);
      for(uint j=0;j<flyCategories.size();++j){
        flyCategories[j]->setMaleAttractionToMales(dNewParam);
      }
      break;
    case 4:
      dNewParam=Flies[0].getMaleAttractionToFemales()+10+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; 
      //for (i=0;i<gciTotalNumberOfFlies;i++)
        //Flies[i].setMaleAttractionToFemales(dNewParam);
      for(uint j=0;j<flyCategories.size();++j){
        flyCategories[j]->setMaleAttractionToFemales(dNewParam);
      }
      break;
    case 5:
      dNewParam=Flies[0].getFemaleAttractionToMales()+10+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; 
      //for (i=0;i<gciTotalNumberOfFlies;i++)
       // Flies[i].setFemaleAttractionToMales(dNewParam);
      for(uint j=0;j<flyCategories.size();++j){
        flyCategories[j]->setFemaleAttractionToMales(dNewParam);
      }
      break;
    case 6:
      dNewParam=Flies[0].getFemaleAttractionToFemales()+10+dChangeAmount;
      if (dNewParam<0)
        dNewParam=0-dNewParam;
      if (dNewParam>dUpperBound)
        dNewParam=dUpperBound-(dNewParam-dUpperBound);
      dNewParam=dNewParam-10; 
      //for (i=0;i<gciTotalNumberOfFlies;i++)
        //Flies[i].setFemaleAttractionToFemales(dNewParam);
      for(uint j=0;j<flyCategories.size();++j){
        flyCategories[j]->setFemaleAttractionToFemales(dNewParam);
      }
      break;
    default:
      throw "\nReach default state of switch in void cExperiment::ProposeNewState(). Exit."; 
  }
  
}

void cExperiment::OutputSummaryStatsAtEndOfSim()
{
    // A list iterator to go backwards in the linked list
    list<cListElement>::reverse_iterator rit;
    // A list iterator to go forwards in the linked list
    list<cListElement>::iterator it;
    //cerr<<"Begin updateSummaryStats\n";
    //BRF_08_20: I added the machinery here to sample the output statistics at equally spaced intervals
    //accounting for the time between moves. I hope.
    int i, j;
    //int  k; 
    float dEndSampleTime;
    //float dStartSampleTime,dSampleInterval, dSampleRemainder;
    //BRF_08_20: These are the variables I need to do the sampling.
    
    // here I'm gonna calculate the actual summary stats (not really 'petri dish' specific,
    // but we'll stick with the lingo - for now)
    //k=0;
    
    for (i=0; i<ciNoStats; i++)
        PetriDishSummaryStats[i] = 0;
    
    // Obtaining end sample time using STL
    dEndSampleTime = statsVector[statsVector.size()-1].dStats[ciNoStats-1];
    //ofile<<"\nendsampletime: "<<dEndSampleTime<<endl;
    
    // ofile<<"\ndEndSampleTime = "<<dEndSampleTime<<"\n";
    // ofile<<"\niSamplingInterval="<<iSamplingInterval;
    
    vector<cListElement>::iterator vit = statsVector.begin();
    // Brad wants to skip the bottom element of the list somehow (BRF_08_20),
    // so we'll do that as well here:
    // rit++; // PM stopped this happening here
    //if (vit!=statsVector.end()){
        //dStartSampleTime = vit->dStats[ciNoStats-1];
        //ofile<<"\nStartSampleTime: "<<dStartSampleTime<<endl;
        //iNSamples-1 intervals gives us iNSamples
        //dSampleInterval=(dEndSampleTime-dStartSampleTime)/(iNSamples-1);
        //dSampleRemainder=0;
    //}
    
    // Now work your way through this list outputting
    vit=statsVector.begin();
    while (vit!=statsVector.end())
    {
        //if ((int)rit->dStats[ciNoStats-1]%iSamplingInterval==0)
        {
            for (j=0; j<ciNoStats; j++)
                ofile<<vit->dStats[j]<<" ";
            ofile<<"\n";
        }
        vit++;
    }
    assert(dEndSampleTime>=vit->dStats[ciNoStats-1]);
}




// PM changed this function
void cExperiment::ReturnToOldState()
{
  // currently we assume all flies have the same parameters - we need to change this function if not
  float eparams[4];
  for(int j=0;j<4;++j){
    eparams[j] = exp(dOldParameterValues[j]/3);
  }
  for (int i=0;i<gciTotalNumberOfFlies;i++)
  {
    Flies[i].RateLeaveMale=dOldParameterValues[0];
    Flies[i].RateJoinMale=dOldParameterValues[1];
    Flies[i].RateLeaveFemale=dOldParameterValues[2];
    Flies[i].RateJoinFemale=dOldParameterValues[3];
    Flies[i].expRateLeaveMale=eparams[0]; //BRF_08_10 // PM added this in case we ever change RateLeaveMale
    Flies[i].expRateJoinMale=eparams[1]; //BRF_08_10 // PM changed these three lines on 3/20/13 - the 3 was missing in the denom
    Flies[i].expRateLeaveFemale=eparams[2]; // PM changed these three lines on 3/20/13 - the 3 was missing in the denom
    Flies[i].expRateJoinFemale=eparams[3]; // PM changed these three lines on 3/20/13 - the 3 was missing in the denom
    //Flies[i].setMaleAttractionToMales(dOldParameterValues[4]); 
    //Flies[i].setMaleAttractionToFemales(dOldParameterValues[5]); 
    //Flies[i].setFemaleAttractionToMales(dOldParameterValues[6]); 
    //Flies[i].setFemaleAttractionToFemales(dOldParameterValues[7]); 
  }
  for(uint j=0;j<flyCategories.size();++j){
    //flyCategories[j]->setMaleAttractionToMales(dOldParameterValues[4]);
    //flyCategories[j]->setMaleAttractionToFemales(dOldParameterValues[5]);
    //flyCategories[j]->setFemaleAttractionToMales(dOldParameterValues[6]);
    //flyCategories[j]->setFemaleAttractionToFemales(dOldParameterValues[7]);
    flyCategories[j]->rollbackParameters();
  }
  if (bToleranceChanged)
    dABCTolerance=gdOldBortotTolerance;
}

// PM added this function
float cExperiment::PriorDensityForTolerance(float dTol)
{
  // Let's use an exponential for now
  return gdPriorLambdaForTolerance*exp(-1*gdPriorLambdaForTolerance*dTol);
}

// PM added this function
float cExperiment::CalculateHastingsRatioForToleranceChange()
{
  //const float cdProbToleranceMovesDown=0.75; // the proble that we propose to substract something from the threshold - note that this could restul in an oncrease to the threshold if we are near the minimum vale for the threshold
  //float gdOldBortotTolerance; // book-keeps the value of the threshold when we are changing it
  //const float gdMaxTolerance=10; // the maximum value the threshold can take
  //const float gdMinTolerance=0.5; // the minimum value the threshold can take
  float dPriorTerm;
  
  // first calculate the ratio of prior terms for the tolerances
  dPriorTerm=PriorDensityForTolerance(dABCTolerance)/PriorDensityForTolerance(gdOldBortotTolerance);
  
  
  // the HR is P(dABCTolerance->gdOldBortotTolerance)/P(gdOldBortotTolerance->dABCTolerance)
  if ((gdOldBortotTolerance-gdMinTolerance)+(dABCTolerance-gdMinTolerance)<gdMaxToleranceChange)
  {
    // we are at the bottom end of the range and a reflection may have occured
    if (gdOldBortotTolerance<dABCTolerance)
      return dPriorTerm*2*cdProbToleranceMovesDown;
    else
      return dPriorTerm*1/2*cdProbToleranceMovesDown;
  }
  else
  {
    if ((gdMaxTolerance-dABCTolerance)+(gdMaxTolerance-gdOldBortotTolerance)<gdMaxTolerance)
    {
      // we are at the top end of the range and a reflection may have occured
      if (gdOldBortotTolerance<dABCTolerance)
      {
        return dPriorTerm*1/2*(1-cdProbToleranceMovesDown);
      }
      else
      {
        return dPriorTerm*2*(1-cdProbToleranceMovesDown);
      }
    }
    else // this change cannot have involved a reflection
    {
      if (dABCTolerance<gdOldBortotTolerance)
        return dPriorTerm*(1-cdProbToleranceMovesDown)/cdProbToleranceMovesDown;
      else
        return dPriorTerm*cdProbToleranceMovesDown/(1-cdProbToleranceMovesDown);
    }
  }  
}

// PM changed this function
void cExperiment::ABC() {
  int iCountOfHowOftenWeHitTheTarget;
  int iNoOfHitsRequired,iNoOfMissesAllowed,iNoOfMissesSoFar;
  float p,dHRterm;
  float dLastAcceptedDistance=500; // will record the distance between observed and simulated data in the last accepted iteration of ABC
  bool bAccept; // records whether we accept or reject each proposal
  // Initialize parameter values
  InitializeEvolutionaryParameters(0); // chose a random set of values - in general we will want to set the argument to 0 // PM did set it to 0 on 1/14/13
  if (bUseBortotScheme)
    dABCTolerance=gdInitialToleranceWhenUsingBortot;
 
  if (mpi_rank==mpi_ranks-1) { 
    ofile<<"\nMPI_Version: Ranseed: "<< ranseed <<"\ntolerance prior param: "
    <<gdPriorLambdaForTolerance<<"  "<<"\nusing Bortot?: "<<bUseBortotScheme
    <<"\nUpperBound for tolerance: "<<gdMaxTolerance<<
    "\nMimunim for tolerance: "<<gdMinTolerance; 
    ofile<<"\nNo of sims per proposal: "<<gciNoOfSimulationsPerProposal<<
    "\nInitial tolerance: "<<gdInitialToleranceWhenUsingBortot;
    ofile<<"\nMale leaving = "<< Flies[1].RateLeaveMale <<"\nMale joining = "
    << Flies[1].RateJoinMale << "\nFemale leaving = "<< 
    Flies[1].RateLeaveFemale << "\nFemale joining = "<< Flies[1].RateJoinFemale;
    ofile<<"\nMale attraction to males = " << 
    Flies[1].getMaleAttractionToMales() << "\nMale attraction to females = " << Flies[1].getMaleAttractionToFemales();
    ofile<<"\nFemale attraction to males = " << 
    Flies[1].getFemaleAttractionToMales() << 
    "\nFemale attraction to females = " << 
    Flies[1].getFemaleAttractionToFemales() << endl;
  }

  if (gciNoOfSimulationsPerProposal%mpi_ranks!=0 || 
  gciNoOfSimulationsPerProposal<mpi_ranks) {
    if (mpi_rank==mpi_ranks-1)
      cerr<<"Error! Would prefer that gciNoOfSimulationsPerProposal\%COMM_WORLD.Get_size()==0  and  gciNoOfSimulationsPerProposal>=MPI::COMM_WORLD.Get_size() for now..."<<endl;
    MPI::Finalize ();
    exit(EXIT_FAILURE);
  }

  float *mpi_dNewSumsSquaresValue=(float*)calloc(mpi_ranks,sizeof(float));
int *mpi_iCountOfHowOftenWeHitTheTarget=(int*)calloc(mpi_ranks,sizeof(int)), *mpi_iNoOfMissesSoFar=(int*)calloc(mpi_ranks,sizeof(int));



  for (int iter=0;iter<gciNoOfIterationsInABCAnalysis;++iter){
    iCountOfHowOftenWeHitTheTarget=0;
    dNewSumsSquaresValue=0;
    // propose new state for fly parameters
    ProposeNewState();
    //cerr<<"proposed new state\n";
    // decide whether of not to accept
    // accept or reject - we accept with prob. iCountOfHowOftenWeHitTheTarget/gciNoOfSimulationsPerProposal
    p=gsl_rng_uniform_pos(rg);
    // PM changed this bit
    // there may now be a term due to the transition kernel for the tolerance
    if (bToleranceChanged) {
      // Calculate the Hastings ratio term
      dHRterm=CalculateHastingsRatioForToleranceChange();
//cerr<<"calculated hastings ratio state\n";
      // We need to multiply our acceptance probability by this. It is equivalent to divide p by it...
      p=p/dHRterm; 
    }
    // we now repeat the simulations only while it is possible to accept
    // How many times do we need to hit the target in order to accept
    iNoOfHitsRequired=1+(int) floor(p*gciNoOfSimulationsPerProposal);
    iNoOfMissesAllowed=gciNoOfSimulationsPerProposal-iNoOfHitsRequired;
    iNoOfMissesSoFar=0;
    // at each iteration we simulate the population gciNoOfSimulationsPerProposal times
    //   j=0;
    // ofile<<"\nIteration "<<i<<"   Number of misses allowed:"<<iNoOfMissesAllowed<<"  No So far:"<<iNoOfMissesSoFar<<" "<<gciNoOfSimulationsPerProposal<<" "<<iNoOfHitsRequired<<" "<<p;
    //  while ((j<gciNoOfSimulationsPerProposal)&&(iNoOfMissesSoFar<=iNoOfMissesAllowed)) {
    //for (j=0;j<gciNoOfSimulationsPerProposal;j++)
    //	j++;

    int iNoOfSimulationsPerRank=gciNoOfSimulationsPerProposal/mpi_ranks;
    //cerr<<"GARY: simulations per rank: "<<iNoOfSimulationsPerRank<<endl;

    for (int k=0;k<iNoOfSimulationsPerRank;k++) {
	if (debugSerialMPI) {
		int global_k = mpi_rank*iNoOfSimulationsPerRank+k;
		long sim_seed = ranseed + global_k; 
		gsl_rng_set(rg,sim_seed);
	}
	// reset values of all statistics
      ResetAllPetriDishStatistics();
      statsVector.clear();
      dSquaredDistance=0, iTimesAFlyMoved = 0, 
      iTimesAMaleMoved = 0, iTimesAFemaleMoved=0;
	// simulate data
      RandomSetup(); // set the flies up randomly
      Evolve(iTimeToStartSampling,0,0); // for this bit we don't record summary stats.
      Evolve(iTotalTimeOfExperiment-iTimeToStartSampling,1,0); // for this bit we do record them, but don't output them

      UpdateSummaryStats();
		 
      if (bUseBortotScheme) {
        if (BortotIsSimulatedDataWithinTolerance())
          mpi_iCountOfHowOftenWeHitTheTarget[mpi_rank]++;//mpi_IsSimulatedDataWithinTolerance[mpi_rank*iNoOfSimulationsPerRank+k]=1;
        else
          mpi_iNoOfMissesSoFar[mpi_rank]++;//mpi_IsSimulatedDataWithinTolerance[mpi_rank*iNoOfSimulationsPerRank+k]=0;
      }else {
        if (IsSimulatedDataWithinTolerance())
	  mpi_iCountOfHowOftenWeHitTheTarget[mpi_rank]++;//mpi_IsSimulatedDataWithinTolerance[mpi_rank*iNoOfSimulationsPerRank+k]=1;
        else
           mpi_iNoOfMissesSoFar[mpi_rank]++;//mpi_IsSimulatedDataWithinTolerance[mpi_rank*iNoOfSimulationsPerRank+k]=0;
      }
    }
	mpi_dNewSumsSquaresValue[mpi_rank]=dNewSumsSquaresValue;

	//put all the random numbers back to the same...
	if (debugSerialMPI)
		gsl_rng_set(rg,ranseed);

	for (int k=0;k<mpi_ranks;k++)  {
		MPI::COMM_WORLD.Bcast(&mpi_dNewSumsSquaresValue[k],1,MPI::FLOAT,k);
		MPI::COMM_WORLD.Bcast(&mpi_iCountOfHowOftenWeHitTheTarget[k],1,MPI::INT,k);
		MPI::COMM_WORLD.Bcast(&mpi_iNoOfMissesSoFar[k],1,MPI::INT,k);
	}

    if (mpi_rank==mpi_ranks-1) {


	iCountOfHowOftenWeHitTheTarget=mpi_iCountOfHowOftenWeHitTheTarget[mpi_rank];
	iNoOfMissesSoFar=mpi_iNoOfMissesSoFar[mpi_rank];

	for (int k=0;k<mpi_ranks-1;k++)  {
		dNewSumsSquaresValue+=mpi_dNewSumsSquaresValue[k];
		iCountOfHowOftenWeHitTheTarget+=mpi_iCountOfHowOftenWeHitTheTarget[k];
		iNoOfMissesSoFar+=mpi_iNoOfMissesSoFar[k];
	}

    if (iNoOfMissesSoFar>iNoOfMissesAllowed)
      bAccept=0;
    else {
      bAccept=1;
      dLastAcceptedDistance=dNewSumsSquaresValue;// not stopping the simulations early anymore  *gciNoOfSimulationsPerProposal/j; // records the average distance, correcting for the fact we may have stopped the simuations early
    }
    
    if (bAccept) {
      // PM changed this bit
      if (!bUseBortotScheme)
      {
        if (dABCTolerance != dABCTolerance){     //this should deal with those bloody NaNs
          dABCNewTolerance = dSquaredDistance;
          dABCTolerance = dABCNewTolerance;
          if (dABCNewTolerance < dThreshold)
            dABCNewTolerance = dThreshold;    
        }else{
          dABCNewTolerance = (dABCTolerance*3 + dSquaredDistance)/4;
          if (dABCNewTolerance < dThreshold)
            dABCNewTolerance = dThreshold;
        }
      }
      iNumberOfAcceptances++;

      // commit parameters now that they're accepted
      for(uint j=0;j<flyCategories.size();++j){
        flyCategories[j]->commitParameters();
      }
      
      // PM add an IF here
      //if (!bUseBortotScheme)
      //  ofile<<"\n  - accepted! The d2dist= " << dSquaredDistance <<" & ";
      
      // PM changed this next bit
      //This bit makes a fixed minimum threshold:
      if ((!bUseBortotScheme))
      {
        if (dABCTolerance <= dThreshold)
          dABCTolerance = dThreshold;
        else
          dABCTolerance = dABCNewTolerance;
      }
    }
    else // reject
    {
      // return to old parameter values
      //  ofile<<"  - rejected ";
      ReturnToOldState();
    }
    
    // PM uncommented this stuff and wrapped it in an if statement, but it probably needs to always apply.
    // output if appropriate
    //if (bUseBortotScheme)
    //{
       if ((iter>gciBurnInPeriod)&&(iter%gciABCSamplingPeriod==0))
      {
        OutputParameterValues(iter,1);
       //cerr<<"GKC DEBUG: last accepted: "<<dLastAcceptedDistance<<endl;
        ofile << " " << dLastAcceptedDistance << "  ";
        OutputSummaryStats();
          ofile<<"     "<<log(Flies[0].expRateLeaveMale)*3<<" "<<log(Flies[0].expRateLeaveFemale)*3<<" "<<log(Flies[0].expRateJoinMale)*3<<" "<<log(Flies[0].expRateJoinFemale)*3;

      }
    //}
  }
}
}

void cExperiment::ResetAllPetriDishStatistics()
{
  iSamplesTaken=0;
  iMoveCounter=0;
  for (int i=0;i<ciNoStats;i++)
    PetriDishSummaryStats[i]=0;
  //I think we need to clear the old petri dishes, too BRF_10_01
  for (int i=0;i<(gciNumberOfPetriDishes+1);i++)
  {
    Dish[i].iHowManyMalesAreHere=0;
    Dish[i].iHowManyFemalesAreHere=0;
  }
}


////////////////THIS IS WHERE PARAMETER VALUES ARE SET

void cExperiment::InitializeEvolutionaryParameters(bool bDeterministic)
{
  // bDeterministic=0; PM removed this on 1/14/13 along with a change to the call of this function from ABC() 
  // if bDeterministic==0 we use random values between 0 and 10; if bDeterministic==1 we use the set of values listed in that part of the the loop below
  int i;
  int iWhichParams=1; // determines which parameters to use when simulating data (and not doing ABC). i.e., when bDeterministic==1
  if (bDeterministic==1)
  {
          
      switch (iWhichParams) // PM added these different settings on 3/1/13 because he wanted to simulated some other datasets
      {
          case 1: // Brad's params
            for (i=0;i<gciTotalNumberOfFlies;i++)
            {                   //0010010
              Flies[i].RateLeaveMale=0;//0.920931;
              Flies[i].RateJoinMale=4;//1.21837;
              Flies[i].RateLeaveFemale=-2;//1.04659;
              Flies[i].RateJoinFemale=5;//1.90516;
              //Flies[i].setMaleAttractionToMales(-6);//0.890104;
              //Flies[i].setMaleAttractionToFemales(5);//2.43081;
              //Flies[i].setFemaleAttractionToMales(1);//1.80603;
              //Flies[i].setFemaleAttractionToFemales(-2);//0.339485;
            }
            for(uint j=0;j<flyCategories.size();++j){
              flyCategories[j]->setMaleAttractionToMales(-6);
              flyCategories[j]->setMaleAttractionToMales(5);
              flyCategories[j]->setMaleAttractionToMales(1);
              flyCategories[j]->setMaleAttractionToMales(-2);
            }
            break;
          case 2: // Brad's params reflect on sex (roughly)
            for (i=0;i<gciTotalNumberOfFlies;i++)
            {                   //0010010
              Flies[i].RateLeaveMale=-2;//0.920931;
              Flies[i].RateJoinMale=5;//1.21837;
              Flies[i].RateLeaveFemale=0;//1.04659;
              Flies[i].RateJoinFemale=4;//1.90516;
              //Flies[i].setMaleAttractionToMales(1);//0.890104;
              //Flies[i].setMaleAttractionToFemales(-2);//2.43081;
              //Flies[i].setFemaleAttractionToMales(-6);//1.80603;
              //Flies[i].setFemaleAttractionToFemales(5);//0.339485; 
            }
            for(uint j=0;j<flyCategories.size();++j){
              flyCategories[j]->setMaleAttractionToMales(1);
              flyCategories[j]->setMaleAttractionToMales(-2);
              flyCategories[j]->setMaleAttractionToMales(-6);
              flyCategories[j]->setMaleAttractionToMales(5);
            }
            break;
          case 3: // Another set that I made up to be different
            for (i=0;i<gciTotalNumberOfFlies;i++)
            {                   //0010010
              Flies[i].RateLeaveMale=2;//0.920931;
              Flies[i].RateJoinMale=2;//1.21837;
              Flies[i].RateLeaveFemale=2;//1.04659;
              Flies[i].RateJoinFemale=2;//1.90516;
              //Flies[i].setMaleAttractionToMales(4);//0.890104;
              //Flies[i].setMaleAttractionToFemales(-3);//2.43081;
              //Flies[i].setFemaleAttractionToMales(-3);//1.80603;
              //Flies[i].setFemaleAttractionToFemales(2);//0.339485;
            }
            for(uint j=0;j<flyCategories.size();++j){
              flyCategories[j]->setMaleAttractionToMales(4);
              flyCategories[j]->setMaleAttractionToMales(-3);
              flyCategories[j]->setMaleAttractionToMales(-3);
              flyCategories[j]->setMaleAttractionToMales(2);
            }
            break;
          case 4: // Those resulting from one of Brad's analyses (ag 6)
            for (i=0;i<gciTotalNumberOfFlies;i++)
            {                   //0010010
              Flies[i].RateLeaveMale=0;//0.920931;
              Flies[i].RateJoinMale=4.05;//1.21837;
              Flies[i].RateLeaveFemale=2.53;//1.04659;
              Flies[i].RateJoinFemale= -2.3;//1.90516;
              //Flies[i].setMaleAttractionToMales(-6.59);//0.890104;
              //Flies[i].setMaleAttractionToFemales(4.33);//2.43081;
              //Flies[i].setFemaleAttractionToMales(1.39);//1.80603;
              //Flies[i].setFemaleAttractionToFemales(1.45);//0.339485;
            }
            for(uint j=0;j<flyCategories.size();++j){
              flyCategories[j]->setMaleAttractionToMales(-6.59);
              flyCategories[j]->setMaleAttractionToMales(4.33);
              flyCategories[j]->setMaleAttractionToMales(1.39);
              flyCategories[j]->setMaleAttractionToMales(1.45);
            }
            break;
      }
      float eRLM = exp(Flies[i].RateLeaveMale/3);
      float eRJM = exp(Flies[i].RateJoinMale/3);
      float eRLF = exp(Flies[i].RateLeaveFemale/3);
      float eRJF = exp(Flies[i].RateJoinFemale/3);
      for (i=0;i<gciTotalNumberOfFlies;i++)
      {                   //0010010
          // PM aded the next 4 lines on 1/23/13 because we were not initializing the exponentiated rate params when simulating data (only)
          Flies[i].expRateLeaveMale=eRLM; // PM divided the exponent by 3 in each of the next 4 lines for consistency with what is below -
          Flies[i].expRateJoinMale=eRJM;
          Flies[i].expRateLeaveFemale=eRLF;
          Flies[i].expRateJoinFemale=eRJF;
      }
  }
  else
  {
    //p=gsl_rng_uniform_pos(rg);

    for (i=0;i<gciTotalNumberOfFlies;i++)  //BRF_08_10
    {
      Flies[i].RateLeaveMale=0;
      Flies[i].expRateLeaveMale=1;
    }
    //Flies[i].RateLeaveMale=2*p;
    //Flies[i].RateLeaveMale=1;
    float probs [4];
    float eprobs [4];
    for( i=0;i<4;++i){
      probs[i] = gsl_rng_uniform_pos(rg);
      eprobs[i] = exp((2*probs[i]/3));
    }
    for (i=0;i<gciTotalNumberOfFlies;i++)
    {
      Flies[i].RateJoinMale=2*probs[0];
      Flies[i].expRateJoinMale=eprobs[0];
      Flies[i].RateLeaveFemale=2*probs[1];
      Flies[i].expRateLeaveFemale=eprobs[1];
      Flies[i].RateJoinFemale=2*probs[2];
      Flies[i].expRateJoinFemale=eprobs[2];
      Flies[i].RateLeaveFemale=2*probs[3];
      Flies[i].expRateLeaveFemale=eprobs[3];
    }
    for(i=0;i<4;++i){
      probs[i] = gsl_rng_uniform_pos(rg);
    }
    for(uint j=0;j<flyCategories.size();++j){
      flyCategories[j]->setMaleAttractionToMales(2*probs[0]);
      flyCategories[j]->setMaleAttractionToFemales(2*probs[1]);
      flyCategories[j]->setFemaleAttractionToMales(2*probs[2]);
      flyCategories[j]->setFemaleAttractionToFemales(2*probs[3]);
    }
  }
  for(uint j=0;j<flyCategories.size();++j){
    flyCategories[j]->commitParameters();
  }
}

// PM changed this function
void cExperiment::OutputParameterValues(int iIteration,bool bOnlyNumbers)
{
  // At the moment all flies are the same, so we just output one of them (Fly 0).
  // One day this will become more complicated
  if (bOnlyNumbers)
    ofile<<"\n"<<iIteration<<"  "<<iMoveCounter<<"  "<<iNumberOfAcceptances<<"   ";
  else
    ofile<<"\nIt: "<<iIteration<<" #jumps: "<<iMoveCounter<<" #accepts: "<<iNumberOfAcceptances<<"   ";

  ofile<<Flies[0].RateLeaveMale<<" ";
  ofile<<Flies[0].RateJoinMale<<" ";
  ofile<<Flies[0].RateLeaveFemale<<" ";
  ofile<<Flies[0].RateJoinFemale<<" ";
  ofile<<Flies[0].getMaleAttractionToMales()<<" ";
  ofile<<Flies[0].getMaleAttractionToFemales()<<" ";
  ofile<<Flies[0].getFemaleAttractionToMales()<<" ";
  ofile<<Flies[0].getFemaleAttractionToFemales()<<" ";
  ofile<< iTimesAFlyMoved<< " " << iTimesAMaleMoved <<" " << iTimesAFemaleMoved <<" ";
  ofile<< dABCTolerance;
}

bool cExperiment::IsSimulatedDataWithinTolerance()
{
  
  bool bTrue=1;
  
  dSquaredDistance = 0;
  //ofile<<"\nSum Stats.: ";
  for (int i=0;i<(ciNoStats-1);i++) //BRF_08_20  Don't want to count the sample times  
  {
    //I changed this to be something like a sum of squared errors
    //dSquaredDistance = dSquaredDistance + pow(PetriDishSummaryStats[i] , 2);
    // GKC:faster to avoid the pow function if exponent is integer;
    dSquaredDistance = dSquaredDistance + PetriDishSummaryStats[i]*PetriDishSummaryStats[i];
    //BRF 09_Mar
  }
  dSquaredDistance = dSquaredDistance + dMAveSD + dFAveSD;
  if (dABCTolerance == 0)
    dABCTolerance = dSquaredDistance;
  dABCNewTolerance = dSquaredDistance;
  
  if (dSquaredDistance > dABCTolerance)
    bTrue = 0;
  
  //cerr<<"GKC DEBUG: squareddist, newssv: "<<dSquaredDistance<<","<<dNewSumsSquaresValue<<endl;
  dNewSumsSquaresValue = dNewSumsSquaresValue + dSquaredDistance/gciNoOfSimulationsPerProposal;
    
    //ofile<<"\ndSquaredDistance= "<<dSquaredDistance<<"  "<<dABCTolerance<<"  "<<dNewSumsSquaresValue<<"  "<<dABCNewTolerance;
  
  return bTrue;
}

bool cExperiment::BortotIsSimulatedDataWithinTolerance()
{
  
  bool bTrue=1;
    
  //ofile<<"\nSum Stats.: ";
  for (int i=0;i<(ciNoStats-1);i++)
  {
    //I changed this to be something like a sum of squared errors
     // dSquaredDistance = dSquaredDistance + pow((PetriDishSummaryStats[i]-TargetPetriDishSummaryStats[i]),2); //PM changed this since we now center the stats
      //dSquaredDistance = dSquaredDistance + pow(PetriDishSummaryStats[i],2);
      // GKC:faster to avoid the pow function if exponent is integer;
      dSquaredDistance = dSquaredDistance + PetriDishSummaryStats[i]*PetriDishSummaryStats[i];
  }
    //BRF 09_Mar
    dSquaredDistance = dSquaredDistance + dMAveSD + dFAveSD;

  if (dSquaredDistance > dABCTolerance)
    bTrue = 0;

  //cerr<<"GKC DEBUG Bortot: squareddist, newssv: "<<dSquaredDistance<<","<<dNewSumsSquaresValue<<endl;
  
  dNewSumsSquaresValue = dNewSumsSquaresValue + dSquaredDistance/gciNoOfSimulationsPerProposal; 
  
  return bTrue;
}

